close all;
clc,clear;
%%  导入数据(时间序列数据一列)
data=xlsread("time_series_data.csv");
data=10*data;
data=data(1:200);
%% 样本
nums = length(data);   %样本数
history = 40;    %历史特征个数

%% 划分数据集
for i=1:nums-history
    re(i,:)=[data(i:i+history-1)',data(i+history)];
end

%%  数据归一化 索引
X=re(:,1:end-1);
Y=re(:,end);
[x,psin]=mapminmax(X',0,1);   %归一化
[y,pcout]=mapminmax(Y',0,1);   %归一化

%% 划分训练集与测试集
num=length(re); %总样本数
% k=input('是否打乱样本(是1，不是0):');
k=1;
if k==0
    state=1:num;
else
    state=randperm(num);
end
r=0.9; %训练集占比
trainnum=floor(num*r);%训练集总样本数
xtrain=x(:,state(1:trainnum));
ytrain=y(state(1:trainnum))';

xtext=x(:,state(trainnum+1:end));
ytext=y(state(trainnum+1:end))';

%% 适应格式
for i=1:trainnum
    x_train{i,1}=xtrain(:,i);
end

for i=1:num-trainnum
    x_text{i,1}=xtext(:,i);
end

%% 创建模型
layers=[
    sequenceInputLayer(history)             %输入层
    lstmLayer(8,'OutputMode','last')        %LSTM层
    reluLayer                               %Relu激活层
    fullyConnectedLayer(1)                  %全连接层
    regressionLayer ];                      %回归层

%% 参数设置
options=trainingOptions('adam' ,...                              %梯度下降算法
                        'MaxEpochs',1200, ...                     %最大迭代次数
                        'InitialLearnRate',0.01, ...             %初始学习率
                        'LearnRateSchedule','piecewise', ...     %学习率下降
                        'LearnRateDropFactor',0.1,...            %学习率下降因子
                        'LearnRateDropPeriod',900, ...           %经过多少次训练后学习率为*0.1
                        'Shuffle','every-epoch', ...             %每次训练打乱数据集
                        'Plots','training-progress', ...         %画出曲线
                        'Verbose',true);                         %观察进度

%% 训练模型
net=trainNetwork(x_train,ytrain,layers,options);

%% 仿真预测
re1=predict(net,x_train);
re2=predict(net,x_text);

%% 反归一化（实际值）
Ytrain=Y(state(1:trainnum));
Ytext=Y(state(trainnum+1:end));

%预测值
pre1=mapminmax('reverse',re1,pcout);
pre2=mapminmax('reverse',re2,pcout);

%% 预测
futurenum=20; %预测数
inputnew{1,:}=y(end-history+1:end);%输入历史数据(进行下一步预测的)
outputs=zeros(1,futurenum); %初始化未来预测值
for i=1:futurenum
%以历史数据作为输入
indata=inputnew{i};
outputs(i)=predict(net,indata(:));
%更新历史数据
inputnew{i+1,:}=[indata(2:end),outputs(i)];
end
%反归一化
reslout=mapminmax('reverse',outputs',pcout);

%% 均方根误差
error1= sqrt(mean(pre1-Ytrain).^2);
error2= sqrt(mean(pre2-Ytext).^2);

%% 相关指标计算
%RE
R1=1-norm(Ytrain-pre1)^2/norm(Ytrain-mean(Ytrain))^2;
R2=1-norm(ytext-pre2)^2/norm(ytext-mean(Ytext))^2;

% MAE
mae1=mean(abs(Ytrain-pre1));
mae2=mean(abs(Ytext-pre2));

%% 图
figure 
plot(1:trainnum,Ytrain,'r-^',1:trainnum,pre1,'b-^');
legend('真实值','预测值');
xlabel('样本点');
ylabel('预测值');
title('训练集预测结果对比');

figure 
plot(1:num-trainnum,Ytext,'r-^',1:num-trainnum,pre2,'b-^');
legend('真实值','预测值');
xlabel('样本点');
ylabel('预测值');
title('测试集预测结果对比');

%百分比误差图
figure 
plot((pre1-Ytrain)./Ytrain,'b-o','LineWidth',1);
legend('百分比误差');
xlabel('样本点');
ylabel('误差');
title('训练集百分比误差曲线');

figure 
plot((pre2-Ytext)./Ytext,'b-o','LineWidth',1);
legend('百分比误差');
xlabel('样本点');
ylabel('误差');
title('测试集百分比误差曲线');

%拟合图
figure
plotregression(Ytrain,pre1,'训练集', ...
    Ytext,pre2,'测试集');
set(gcf,'Toolbar','figure');

%未来预测图
figure
plot(1:nums,data,'r-o','LineWidth',1);
hold on
plot(nums:nums+futurenum,[data(end);reslout],'b-*','LineWidth',1);
xlabel('时间');
ylabel('值');
legend('历史数据','未来预测数据');
grid on














